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Abstract 

As a departure from existing continuum approaches for describing the stabihty 
and evolution of surfaces of crystalhne materials, this article provides a description 
of surface evolution based on the physics of the main feature imposed by the discrete 
nature of the material, namely, crystallographic surface steps. It is shown that the 
formation energy of surface steps depends on the sign of extensional strain of the 
crystal surface, and this behavior plays a crucial role in surface evolution. The nature 
of this dependence implies that there is no energetic barrier to nucleation of islands 
on the growth surface during deposition, and that island faces tend toward natural 
orientations which have no counterpart in unstrained materials. This behavior is 
expressed in terms of a small number of parameters that can be estimated through 
atomistic analysis of stepped surfaces. The continuum framework developed is then 
applied to study the time evolution of surface shape of an epitaxial film being deposited 
onto a substrate. The kinetic equation for mass transport is enforced in a weak form 
by means of a variational formulation. It is found that islands form without nucleation 
barriers and they evolve to shapes with natural surface orientations. The implications 
of the calculations are shown to be consistent with the behavior observed during 
deposition of semiconductor materials in recently reported experiments. Finally, it is 
verified that the predictions of the continuum model are essentially the same as those 
of the discrete step model for an isolated strained island. The development in this 
article is limited to two-dimensional plane strain deformation to keep the arguments 
transparent, but this is not a fundamental limitation of the approach. 

Keywords: surface diffusion, surface energy, morphology evolution, semiconductor 
material, stability and bifurcation 

1 Introduction 

Self-organized semiconductor nanostructures hold the promise for manufacture of micro-electronic 
devices with unprecedented performance characteristics. Because of the close connection between 
size and electronic characteristics, these structures can be tuned to very specific requirements. 
Potential applications include field effect transistors, quantum memory devices, and solid-state 
lasers. 



Strain-driven nucleation, growth and coarsening of epitaxial islands offers a versatile ap- 
proach to manufacturing nanoscale devices. To exploit this phenomenon in manufacturing, a 
fundamental understanding of the role of strain, surface energies, and the kinetics of transport 
on the formation and evolution of material structures is required. 

Surface energy is an important concept in considering the evolution of microstructure 
in small-scale material systems. The free energy of the bounding surface of a crystal is a 
macroscopic quantity representing, in some sense, the net work that had to be done to create 
that surface. This energy is distributed over a specific mathematical surface that approximates 
the physical boundary between the crystal and its surroundings. Prom this definition of surface 
energy it is clear that the reference level for energy of a free surface is the state of the material 
on that same crystallographic surface when it is embedded deep within a perfect crystal. A 
number of analyses of surface stability and evolution have been reported. Asaro and Tiller 
(1972), Grinfeld (1986) and Srolovitz (1989) independently showed that the flat surface of a 
stressed solid under two-dimensional plane strain conditions is unstable if the surface energy is 
independent of surface orientation. The same physical model was extended to three dimensions, 
but restricted to small amplitude surface fluctuations, by Freund (1995). A numerical method 
for handling large amplitude surface fluctuations was developed by Zhang and Bower (1998). 
Gao and Nix (1999) adopted the same physical framework for describing the breakup of an 
unstable film into islands. In all of these studies, the surface energy was assumed to be isotropic, 
independent of surface strain and to have no connection to the discrete nature of the crystalline 
materials being modeled. This study is focused on a departure from this point of view. 

The value of surface energy per unit area of a given crystallographic surface orientation 
is determined by the fine scale structure of that surface. For a high symmetry orientation in 
a crystal, such as a {100} surface of a cubic crystal, the surface is atomically flat. For other 
orientations close to this surface, the structure usually consists of flat terraces with well defined 
local surface energies, separated by atomic scale ledges or steps as illustrated schematically in 
Fig. ||. The steps alter the macroscopic surface energy by an amount corresponding to their 
energy of formation in the configuration relevant to the structure. Furthermore, the energies 
of individual features and the interactions of these features are influenced by the presence of 
strain in the crystal. Macroscopically, the surface is assumed to be smooth but to have a local 
surface energy at a point on the surface determined by the orientation of the tangent plane at 
that point and the level of elastic strain in the crystal. A quantitative interpretation of the 
free energy of a strained crystal in terms of the physics of crystallographic steps is developed in 
the sections that follow. The driving force for alterations in surface morphology is the surface 
chemical potential, defined in terms of variations of system free energy with surface shape. The 
surface chemical potential of a stepped surface is also derived, and implications for morphology 
evolution by surface diffusion are examined. The actual values of the parameters involved in 
the characterization can be determined only through atomistic simulation or experiment; such 
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work focused on the SiGe material systems is being carried out in parallel with this study. 




Figure 1: Geometry of a vicinal surface that makes an angle of 9 with a high 
symmetry direction. The spacing d between the steps is related to the step-height hg 
through the relation d = hg/ taxiO. The dashed line shows the macroscopic surface 
orientation. 



The present work is motivated by recent experimental observations reported by Sutter and 
Lagally (2000) and Tromp et al. (2000) on the growth of Sii-^Gea; films on Si(OOl) substrates 
with Ge concentration x in the range 0.1 < x < 0.4. The experiments have clearly revealed 
that, during early stages of growth, shallow stepped mounds whose side walls are made up 
of widely spaced steps emerge as an inherent morphological instability of the film. As more 
material is deposited, the spacing between the steps gradually decreases until the sidewalls reach 
a certain crystallographic orientation. This type of nucleationless growth of epitaxial islands is 
not predicted by a widely used nucleation barrier model for the energetics of faceted islands 
(Tersoff and LeGoues 1994), which is based on the competition between gain in elastic energy 
and the energy spent in creating the surface of the island. Approaches to surface instability that 
account for neither orientation dependence nor strain dependance of surface energy have been 
unable to explain the observed growth mode. A similar growth mode has also been observed in 
very recent experiments on Ge films grown on Si(OOl) substrates (Vailionis et al. 2000, Rastelli 
et al. 2001), where small islands with widely spaced steps called " pre- pyramids" are the three 
dimensional features that first appear during growth; these pre-pyramids eventually evolve to 
faceted islands. In this paper, we will show that the emergence of stepped islands, without a 
nucleation barrier, can be explained by incorporating the physics of steps, in particular their 
interactions and the dependance of their formation energies on the mismatch strain. 

Before we proceed with development of the continuum description, we recall a few results 
regarding stepped crystal surfaces that are well-established in surface science. Below the char- 
acteristic roughening temperature, which is in excess of 1000 °C for most semiconductors, a 
nominally fiat surface of a crystal that is misoriented by a small angle from a high-symmetry 
direction consists of a train of straight parallel steps, as illustrated schematically in Fig. |^. The 
surface energy of such a "vicinal" surface is given by 

I sin 

7(0) =7ocos0 + /3i|sin0| 2^, (1) 
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where 6 is the misorientation angle, 70 the energy density of the atomicahy flat surface and 
the parameters /3i and P3 are related to step creation and interaction energies respectively. 
Specifically, /3i = l3/hs, where /? is the energy to create a unit length of an isolated step and hg 
is the height of an atomic step. Like any other defect, say a dislocation or a vacancy, steps give 
rise to long-range stress fields in the solid; the steps interact with each other, as well as with 
other defects in the material, through these fields. The amplitude of the stress field produced 
by a step decays with the inverse square of the distance from the step (Marchenko and Parshin 
1980). The last term on the right in ([l|) arises from the interactions between steps through their 
stress fields. 

The discussion proceeds in the following way. In Sec. 2 we will derive the energy of a 
stepped surface for small strains but arbitrary orientations, starting from the surface energy of 
a vicinal surface given in (^). To facilitate calculations, this energy will then be simplified so 
that it applies to surfaces whose deviation from a flat orientation is small. The surface chemical 
potential which represents the driving force for change in shape due to mass transport on the 
surface will be derived in Sec. 3. In Sec. 4, we develop a variational formulation of surface 
evolution and apply it to study the growth of a thin film on a lattice mismatched substrate 
under the influence of a constant deposition flux. Our focus will be on the evolution and 
interactions of stepped mounds that appear on the substrate surface as soon the deposition 
flux is turned on. It will be shown that the spacing between the steps on the sidewalls of 
these mounds decreases continuously until it reaches an optimum value, determined by the 
competition between repulsive step interactions and a strain induced lowering of step formation 
energies. We provide a summary of the key results along with future directions of research in 
Sec. 5. The connections between the continuum description developed in this paper and discrete 
elastic models of steps will be discussed in Appendix A. In Appendix B, we calculate the energy 
of a strained island using both the continuum and discrete descriptions of stepped surfaces. 
Here, we show that the mismatch strain can lead to lowering of the surface energy of the vicinal 
surfaces and can lead to elimination of the barrier to nucleation of these islands. 

2 Energy of a strained solid with a stepped surface 

Consider a strained crystal that occupies the region R and that undergoes two-dimensional 
generalized plane strain deformation. The region is bounded by a free surface S as shown in 
Fig. |2|; the remainder of the boundary of R is "workless" due to symmetry constraints or its 
remoteness. The free surface is viewed as a being made up of an array of inflnitely long, straight 
monatomic steps of height hg . Initially, the surface S is atomically flat with its normal vector in 
the X2— direction. The goal of this section is to represent the change in energy associated with 
departure from flatness in terms of the evolving surface shape. 

When the surface S is flat, the crystal is subjected to a spatially uniform strain e^*-. To 
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ensure that S is stress-free in this configuration, this initial strain must satisfy 

Ci2ki<^ki = 0' (2) 

where Cijki is the array of elastic constants. The total strain in the crystal is then written as 
^ij ~ "I" + ^ij' that Ui represents the additional displacement field that arises as 

the surface deviates from its flat orientation. The material stretch along the evolving surface 
can be expressed in terms of the limiting value of the bulk strain on the surface as e = mieijUij, 
where rui is a unit vector that is locally tangent to S. For the surface of the stressed solid under 
consideration, as illustrated schematically in Fig. |^ we can generalize (|l]) to include the effects 
of surface stretch e to express the surface energy density as 

7[h{s),e{s)] = {jo + roe{s))^l-h'{sr + {Pi + Pie{s)) \h' {s)\ + (3,,^^^^, (3) 

where h{s) is the deviation of the height of the surface from its fiat orientation as a function of 
arclength s along S, tq = 70 + djo/de is the surface stress of the flat surface at the current level 
of stretch, and (3i = df3i/de is a measure of the sensitivity of the formation energy of a step on a 
strained surface. The last term has the same interpretation as in (||) and its inclusion is essential. 
In principle, we should also expand the coefficient of |/i'(s)|^ in terms of surface stretch. However, 
it will be shown in the appendix, on the basis of direct comparison with discrete step models, 
that ignoring the higher order terms in the expansion provides a reasonable approximation. 
Consequently, only the leading order term represented by j3^ is retained. 




Figure 2: Configuration of the stepped surface 5 of a strained crystal R. The 
deviation of the height of the surface from its initially fiat orientation is denoted by 
h[s) where s is the arclength of the surface along S. 

The total energy of the solid is the sum of the elastic energy and the surface energy, 

^[/i(s),eij(x)] = \ I Cijkieij{yi)eki{yi)dy. + I Y^{s)ds + f T{s)e{s)ds, (4) 
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S(.) = 70^1 + A \h'{s)\ + /33 Y37^ (5) 

and 



where x is any material point in R, and and t(s) are the local surface energy and surface 
stress at zero surface stretch. The latter quantities are given by 

\h'{s)l 

T{s)=ToyJl-h'{s)^+/3i\h'{s)\, (6) 

respectively. Because mechanical equilibrium is achieved very quickly on the time scale as- 
sociated with diffusive mass transport on crystal surfaces, we will first derive conditions for 
equilibrium at fixed surface shape h{s). Then, by using the equilibrium fields obtained, we can 
proceed to derive an expression for the change in the energy of the system in terms of surface 
shape change, which leads naturally to a corresponding definition of surface chemical potential 
as an energetic driving force for surface shape change. 

Mechanical equilibrium rcqiiires that the variation in the free energy E at fixed surface 
shape h{s) due to a small perturbation in displacement field 5ui from equilibrium values of Ui 
must vanish to linear order in 6ui. The first order variation in energy in the present case is 

6E = - Cijkiui^kj^Uidx. + / Cijki (ui^k + ^kl) rijSuids + / T{s)mi6uijmjds, (7) 

J R J S J S 

where nj is the local outward unit normal vector to S. The last term can be conveniently 
rewritten as 



i 



f f \ dSui 

T{s)miduijmjas = J T[s)mi— — as. (8) 



IS 

Because the bounding surface of R is workless, we can rewrite this expression in the form 

/ ris)m,^ds = - [ ^^^^^6u,ds. (9) 
Js ds Js ds 

Noting that 

d(T(s)mi) dris) / ^ / x .-^^ 

^ V ' = —r^rrii + T ni, 10 

ds ds 

where k{s) is the curvature of the surface, which is taken as positive if the center of curvature 
is outside the material, we can rewrite the variation in the total energy as 

dT{s) 



SE = - Cijkiui^kj^uidyi + Cijki {ui,k + nj ^^"^^ " T{s)K{s)ni 



Suids. (11) 



Because the variation must vanish for arbitrary 6ui, we obtain the usual equilibrium condition 
that the divergence of stress must vanish at each material point, CijkiUi^kj = in the region R. 
In addition, traction boundary conditions reflecting the structure of the surface emerge in the 
form 

Cijkm^knjrrii = -Cijkielinjrrii + /3i/t"(s)Sgn[/i'(s)] - ^^^if^=^, 

VI - h'{sY 

Cijkm,knjni = -CijkielinjHi + (^tq^Ji - h'{s)'^ + 0i \h'{s)\j k{s), (12) 
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where Sgn[p] = d\p\/dp for any quantity p 7^ 0. These must be satisfied on S by the elastic 
field. The first terms on the right sides of (12) represent, respectively, the normal and shear 
tractions applied to the surface to compensate for the relaxation of the bulk stress Cijkie^j due 
to surface reorientation; these tractions ensure that the surface S remains traction free. The 
other terms on the right side of ([T2|)i represent the shear traction due to variation of surface 
stress with position along the surface, whereas the second term in (|l2|)2 represents the normal 
traction on S due to the so-called Laplace pressure which arises from surface stress. Note that 
these conditions arise only as a result of the assumption that energies of surface features may 
depend on strain. They involve no other ad hoc assumptions. 

Up to this point, there has been no restriction on the magnitude of the surface shape 
change. From this point onward, however, attention will be restricted to surface shapes that 
are locally misoriented by only a small angle from the flat orientation, that is, the condition 
<C 1 prevails everywhere on S. Using the fact that mi = 1 — h'{xi)'^/2 = n2 and 
771-2 = h'[xi) = —rii to second-order in the slope, the boundary conditions become 

Ci2kiui,k = Ciikielih'{xi) + 0ih"{xi)Sgn[h\xi)] = fi{xi), 

C22klUl,k = TQh"{xi) = /2(xi), (13) 

where /i and /2 are the components of surface traction in the coordinate directions that give 
rise to stress fields in the crystal associated with shape change. The displacement field at any 
point in the crystal can now be obtained using the elastic half-space Green's function Gij{x) as 

/oo 
Gijix- x'ei)fj{x')dx'. (14) 
~oo 

We can now recast the minimum free energy at fixed surface shape in terms of the surface 
tractions introduced in (13). Applying the divergence theorem to the first term in (Q), recalling 
the workless nature of the boundary constraints, and using the fact that 

(to + Pi m,e° mj = (tq + (3i eJi + 2Toe^^2h' 

+ (^o(e^2 - 4) + 2/3ie;2Sgn [h']) {hf + 0{h^) (15) 

we find that 

E[h{x)] = Eo + J^ (/3i+/3ie?i + 2Toe?2Sgn[/i'(x)]) |/7'(x)|(ix + ^ /?3|/i'(x)|^dx 

+ (to(4 - en) + 2/3ie?2Sgn [h'ix)]) {h'{x)fdx - ^ fi{x)u,{x)dx, (16) 

where Eq is the total energy when the surface is fiat, Si refers to the projection of S in the 
xi-direction and the coordinate xi has been replaced by x for simplicity. This is the total free 
energy of the system with a perturbed surface shape. 

Before proceeding to study the evolution of stepped surfaces, we make a few observations 
about the continuum representation of system free energy change given by (^) and (16): 
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1. All the terms in ( |16D except the last are local, in the sense that the energy density at 
any point depends only on the local surface slope. The last term, however, represents the 
interaction of the surface force distributions given in (|l^); in the jargon of surface science, 
these represent force monopoles. Since these monopole interactions are long-ranged, the 
energy density at a point on the surface depends on the surface shape at all points on the 
surface. 

2. The free energy in the small slope approximation is correct to 0{h?); only one term of 
0{h^) has been retained. It will be shown in the Appendix A that this approximation 
where all the other 0{h^) terms are ignored is a reasonable one, by means of comparison 
of the continuum description with discrete step models. 

3. The small slope approximation is not essential to the development. If this approximation is 
not made, the total energy can be evaluated by solving the boundary value problem, which 
consists of the bulk equilibrium equations subject to boundary conditions (12), using a 
numerical technique such as the finite element or boundary element method. 

4. The quantity f3i must be either positive or negative, if the step formation energy is strain 
sensitive. It follows that /3i + /Sie^i < for compressive or tensile strain e^i, respectively. 
This implies that the energy of the stepped surface becomes lower than the energy of the 
surface without steps for strain of one sign or the other. A plot of the surface energy of a 
vicinal surface for which = is shown in Fig. |3|. In this case, when strains are sufficiently 
compressive with [3i > 0, the variation of surface energy with surface orientation 0, as given 
by 

liO) = 70 + roe?i + + Piel,) \e\ + roie^^ - e?i)^' + Al^l', (17) 

develops minima away from 9 = 0. At small misorientation angles, the step formation 
energy dominates the repulsive interaction energy between the steps, while at larger angles, 
the step interactions are larger in magnitude. The competition between these two opposing 
effects results in optimum misorientation angles denoted by 9 = ±0* in Fig. ^. It follows 
that Pi is a key quantity in determining the morphology of strained crystals^. Interestingly, 
since the force monopole in (|^) has a contribution that depends on this key quantity 
also enters the non-local part of the total energy. 

5. If the energy of a stepped surface under strain becomes lower than the energy of a flat 
surface, there is no barrier to nucleation of epitaxial islands whose sidewalls are made 
up of surfaces with widely spaced steps. In Appendix B, we consider the energetics of 



^The value of /3i cannot be deduced from a continuum analysis. However, by means of atomistic simulations 
of stepped surfaces on strained crystals of Si being conducted in parallel with this study, we have concluded that 
the anticipated minima do indeed arise within the strain range of interest. 
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such stepped mounds using both the discrete models of steps and the present continuum 
description. 



Tensile strain 
, Unstrained 




-e* 9* 



Figure 3: Dependence of the surface energy 7(0) of a vicinal surface on surface ori- 
entation Q for /?! > 0, To > and ^2 — 0- When compressive strains are sufficiently 
large, the energy of a stepped surface becomes lower than the energy of a flat surface. 
In this case, competition between strain-induced lowering of the formation energy 
of the steps and the repulsive step interactions results in a minimum in the surface 
energy at an optimum angle, Q* indicated in the figure. 



3 Surface chemical potential 

The evolution of surface morphology of strained crystals as a result of mass transport by surface 
diffusion or other mechanism has been studied extensively on the basis of continuum models, 
beginning with the introduction of a surface chemical potential by Herring (1953). The essential 
ingredients for describing evolution are the chemical potential field defined as a function of 
position over the evolving surface of the solid and a kinetic relationship between surface mass 
flux and the gradient in chemical potential. The rate of change of surface shape then follows from 
conservation of mass. The kinetic relationship must be defined in such a way as to ensure that 
the free energy decreases as surface evolution progresses. In this section, the surface chemical 
potential for a stepped surface will be obtained from the free energy expression given in (^. 
The result is written first for arbitrary magnitude of surface slope, and it is then reduced to a 
simpler form applicable for the case of surface perturbations with slopes of small magnitude. 

In deriving the mechanical equilibrium equations in Section 2, it was assumed that the 
surface was fixed with respect to the material. This was justified by the observation that 
mechanical fields equilibrate very fast on the timescale of surface evolution. We now assume 
that the surface shape depends on time t, and we consider the associated rate of change the 
free energy of the system as the surface shape changes, with the system being continuously in 
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mechanical equilibrium. The rate of change of this free energy measure then has the form 

E{t) = J^fi[h{s,t)]vn{s,t)ds, (18) 

where Vn{s,t) is the outward normal velocity of the surface S with respect to the material 
instantaneously on it and n[h{s, t)] is the local chemical potential. At any point along the 
surface, the quantity Vn ds is the rate of addition of material volume to the surface. It follows 
that its coefficient is the change in free energy per unit volume of material added, that is, the 
coefficient /U represents the surface chemical potential field (Freund 1998). 

We can write the rate of change of free energy (|^ , following application of the divergence 
theorem to the first term, as 

^ f rr . r [ f 8^ dO Bt dO cle \ , [ ,^ 

Js ^ Jg yoe'dt 'dO'dt^^ ) Js KVnds, (19) 

where U = ^Cijkieijeki is the elastic strain energy density, a superposed dot denotes a material 
time derivative, and d /dt denotes the total time derivative whenever it must be distinguished 
from the material derivative. For example, the rate of change of surface stretch 



de 
dt 



— = e + rmmjeij^knkVn (20) 



includes a convective term, and the rate of change of surface orientation following the non- 
material surface requires that 

d6 dvn , drrii dvn 
11 = ^ ^ = "^^- 

The first term in (^) accounts for the change in strain energy in the crystal due to surface 
evolution, whereas the remaining terms account for the change in surface energy. Recalling that 



the system is enclosed within workless boundaries, we can rewrite (1£) in the form ([Lq) with 

/i[/l(s, t)] = U — K 



dr de d . , , , 

d6d^~ ^'ds^'^''^^^'^^' + mimj-ey- fcUfcT. (22) 



In the limit of small slope \h^x (^^lOI 1> where the subscript x following the comma indicates 
partial differentiation with respect to x, this expression can be further simplified to the form 

/x[/i(x, t)\ = Uo- h,xx [(/3i + pie''u)6 {h,, ) + 2ro(e^2 " e?i) + ^Pi^u^gn [h,^ ] + 6/^3 \h,, \ 

+CukieliUi^x - ASgn [h (23) 

where Uq is the value of U when the surface is flat and 6{ ) is the Dirac delta function. 



A vexing feature of the chemical potential function obtained in (23) is the appearance of 
the Dirac delta function in the term representing step creation energy. We expect this term to 
play a role in evolution of the surface of the strained crystal because the step creation energy can 
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become negative for strain of one sign or the other, if the magnitude is sufficiently large. While 
it is possible to introduce any number of ad hoc "regularizations" to deal with this singularity, 
it is not evident that the results obtained will be independent of the particular choice made. 
Therefore, we adopt a variational framework for describing surface evolution which is based on 
the chemical potential derived but which circumvents the problem posed by the singular step 
creation energy term. 

4 A variational approach for strain-driven surface evolution 



In this section, we apply the continuum description of the energy derived in (16) to study 
the growth of a strained film bonded to a lattice-mismatched substrate. The constraint of 
the substrate maintains the mismatch stress in the film. As noted before, recent experiments 
(Sutter and Lagally 2000 and Tromp et al. 2000) have observed that growth of the deposited 
film proceeds via formation of shallow mounds whose stepped side-walls ultimately evolve to a 
faceted orientation. The evolution takes place through diffusion of atoms across the terraces 
that separate the steps on the side- walls of these mounds. 

In the case of terrace diffusion-limited kinetics, the mass flux on the surface can be related 
to the gradient of the chemical potential derived in (23) through 



j{x,t) = -Cfl^^r, (24) 

where c is a coefficient representing surface mobility of diffusing species. If there is a constant 
deposition flux /, mass conservation can be invoked to connect the mass flux to the surface 
shape through the relation 

ht+j,x = f. (25) 

Our goal is to develop a variational framework for modeling the evolution of these stepped 
surfaces. If we focus attention on surface profiles that are periodic with wavelength A, the 
thickness of the film h{x,t), which is illustrated in Fig. ^, can be expressed as 

Nf 

h{x, t) = ^ an{t) cos{nkx), < x < \, (26) 

n=0 

where k = 27r/A. The wavelength should be large enough to accommodate several stepped 
mounds and other features of interest. The number of terms Nf included in the series is finite 
but otherwise unrestricted. The variational formulation provides first order coupled ordinary 
differential equations for the Fourier coefficients a„(t); these equations will be numerically inte- 
grated to obtain the evolving surface shape for t > using information on the surface shape at 
t = 0. 
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Using (Pq) in (|25|), we can express the mass flux in the form 

sin(nA;x) 



Nf 



j{x,t) = - ^ dn{t)- 



n=l 



nk 



(27) 



and the mean height of the surface as ao(^) = ft- While (24) and (25) constitute a closed set of 
equations, they cannot be directly integrated to obtain the surface shape because of the singular 
nature of the chemical potential in (|2^), as noted in the preceding section. We now turn to the 
variational approach which provides well-behaved evolution equations for the Fourier coefficients 
in (p^) by using the weak form of (|23). 





h(x,t) , 






Film 


r ' 



Substrate 



Figure 4: Strained film on a relatively thick lattice-mismatched substrate. The time 
dependent thickness of the film is denoted by h{x,t). When the film surface is fiat, 
there is a spatially uniform elastic mismatch stress in the film and the substrate is 
stress free. 



The variational approach can be formulated (Suo 1997) by noting that the functional 

^j{x,t)] = [\,,j{x,t)dx+ f^Li^dx (28) 
Jo Jo 2c 



attains a minimum value when j{x,t) satisfies the kinetic relation (24). The idea is to express 
the functional in terms of the time derivatives of the Fourier coefficients d„(t) and to then 
minimize the functional with respect to variations in these rates. The first of these tasks can be 



accomplished by integrating the first term by parts and by using (|25|) and (^) to obtain 

Nf Nf 
n=l " n=l 

where the total free energy E[a(),ai ■ ■ ■] of the film-substrate system can be directly evaluated 
from (16) with one minor modification. Since the mismatch strain is present only in the film 
but not in the substrate (refer to Fig. ^), e^j in ( p^ ) has to replaced by a spatially dependent 
function e^j{x2) = e^jH{x2), where H{x2) is the Heaviside function which has the value of zero 



12 



when rE2 < and the value one when X2 > 0. Minimizing the functional in (|29| ) with respect to 
the d„'s, we get the evolution equations 

cnk^ dE 

an{t) = T^, 30 

TT dan 

which are unambiguous and well-behaved. We now proceed to integrate these rate equations 
to determine the time evolution of the deposited film. At each time step of integration, the 
quantity dE/dan is calculated numerically. The important role of the nonlinear contribution 
to ( |30[ ) for the convergence of this approach is illustrated in Appendix C by considering only a 
single Fourier mode. 

The numerical integration of ( |30| ) was carried out using the fourth-order Runge-Kutta 
procedure with adaptive step-size control as described in detail by Press et al. (1992). The 
parameters in the continuum description of the energy were chosen as follows: /3i = 0.03 J/m^, 
/?! = 15J/m?, tq = IJ/m^, e^i = —0.01 and (5^ = 2.8QJ/m?. It is assumed that the substrate 
and the film are isotropic with similar elastic properties, and that the Young's modulus and 
the Poisson's ratio are 10^^ N/m? and 0.3, respectively. Using these parameters, we find that 
^22 ~ 0.0023 and = 0. Also, for the compressively strained film, the surface energy of the 
film (sketched schematically in Fig. ^ attains a minimum when 9* = 0.12, which implies that 
the sidewalls of the stepped mounds would eventually evolve to this angle. 

The initial profile of the film was chosen to be a sinusoid of wavelength A = 400nm, with 
an amplitude of 0.4 nm, so that the only non-vanishing Fourier components at t = are oq 
and ai. In our calculations, sixteen Fourier coefficients were used to keep track of the surface 
shape. The evolution of the deposited film is shown in Fig. ^. As the deposition flux is turned 
on, the material on the surface very quickly gathers into five stepped mounds with slopes much 
smaller than the optimum value of 0.12 as shown in Fig. ^(a). The insert in Fig. ^(a-c) tracks 
the evolution of one of these mounds. As more material is deposited, the side-walls of the mound 
become steeper until they reach the optimum angle ^ = 0.12 in Fig. |5|(c). It can also be observed 
that the center of mass of the mound shifts gradually to the right in going from Fig. |5|(a) to 
Fig. ||(c). This can be understood by looking at the interactions between the islands. Since 
elastic relaxation is achieved for widely spaced islands, the islands tends to repel each other 
(Floro et al. 2000). It can be seen that the island in Fig. |5|(a) is located closer to the island 
on its left and would therefore tend to shift towards the right via diffusion of atoms from the 
side- wall on the left to the one on the right. Once the side- walls of the islands reach the optimum 
orientation, they grow in a self-similar fashion, until they come in contact with their neighbors 
as shown in Fig. |5|(c). At this point, self-similar coarsening is initiated, which leads to a decrease 
in areal coverage (fraction of the surface of the substrate covered by the film) as is evident in 
Fig. ^(e). Here there are two islands, with side walls a,t 9 = 0.12, separated by 50nm. As more 
material is deposited, these further grow in size and, eventually, one of the islands grows at the 
expense of the other island as shown in Fig. ^(f)-(h). 
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Figure 5: The figure shows a time sequence of surface profiles of h{x, t) plotted on the 
vertical axes versus x plotted on the horizontal axes of a strained film during constant 
flux growth on a lattice-mismatched substrate. All the dimensions in nanometers. 
The inserts in (a)-(c) show the evolution of the third island from the right. To aid 
in the comparison of shapes at different times, we have included the island shape 
from (a) in (b). Similarly, we have included the island shapes from (a) and (b) in 
(c). The slope of the largest island in each of the smaller insets is indicated in the 
figure. 
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The evolution of the film in the early stages of island growth are in close agreement with 
the recent experiments of Sutter and Lagally (2000) and Tromp et al. (2000). The key result 
of these experiments is that the islands evolve as a natural instability without any nucleation 
barrier. As we show in the Appendix, there is no barrier for nucleating stepped islands if the 
strains are compressive. Islands with orientations below the optimum facet orientation lower 
the energy of the film and provide a kinetic pathway for obtaining faceted islands that is free 
of any nucleation barrier. This is indeed what is seen in experiments and during early stages of 
growth in Fig. ^(a-c). 

It can be seen from Appendix B that the energy of the islands is dominated by the 
surface energy at small island volumes and by the elastic energy at large volumes. Using the 
parameters adopted in our calculations, we can use (43) to show that the crossover between 
these two regimes takes place when the base width of the island is about 200-300nm. When the 
island sizes are smaller than this value, the sidewalls are oriented at the optimum angle that 
minimizes the surface energy. With increasing base width, it is known that the islands undergo 
shape transitions (Medeiros-Ribeiro et al. 1998), where the sidewall angles change to a steeper 
orientation, which is usually a low-energy crystallographic orientation that makes a larger angle 
with the substrate. Since such orientations are not yet included in our model, we do not observe 
such transitions. We do, however, find that the slope of the sidewalls of the large islands in 
Fig. |^(e) is about 0.15, which is about 25% bigger than the optimum slope. This indicates that 
elastic energy of these larger islands is becoming comparable to the surface energy. 

There is no fundamental impediment to including additional low energy surface orienta- 
tions in order to study steeper sidewall facets. For example, these might appear as additional 
relative minima in the variation of ^{0) with orientation 9 at angles greater in magnitude than 
9*. This will be subject of future investigations. 



5 Summary 

In summary, we have developed a continuum framework to model the energetics and evolution 
of stepped surfaces of nanostructures. The surface energy of strained surfaces is obtained by a 
generalization of the energy of a vicinal surface to account for the effects of stretching of the 
surface caused by the mismatch strain. We find that the dependence of the formation energy 
of steps on the sign of mismatch strain plays a crucial role in establishing the morphology of 
evolving nanostructures. While this dependence has been seen in atomic-scale simulations, we 
have developed a consistent continuum framework that naturally accounts for this crucial effect 
in modeling surface evolution. Furthermore, we have shown in Appendix A that our continuum 
description provides a direct way to analyze stepped surfaces, without evaluating complicated 
sums involving discrete steps. 

The continuum framework was applied to study the evolution of an epitaxial film bonded 
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to a lattice-mismatched substrate. Our simulations show that the deposited material initially 
gathers in shallow stepped mounds whose side- wall angles eventually evolve to an orientation that 
is determined by the competition between strain-induced lowering of the step formation energy 
and the repulsive interactions between the steps. This kinetic pathway has no nucleation barrier 
and occurs as a natural instability, in agreement with the recent experimental observations. The 
simulations were also used to look at self-similar growth of islands leading to their impingement 
with neighboring islands and subsequent coarsening. 

In this article, discussions were limited to two-dimensional deformation fields and evolu- 
tion of surfaces in one dimension through the diffusion atoms, all of the same type. Work is 
currently in progress to extend the current formulation to study the evolution of two-dimensional 
surfaces. Also, since the films of interest in device applications are typically alloys of two dif- 
ferent semiconductors, further insight into surface evolution can be obtained by including the 
possibility of alloy segregation in binary materials. 
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Appendices 



A Discrete surface step models 

In developing a continuum description of stepped surfaces, we started with the surface energy 
of a strained vicinal surface given by (^, which is a generalization of the surface energy of 
the vicinal surface to include the effects of surface stretch. The key idea of this section is to 
demonstrate that if one chooses to start with elastic models of discrete steps and proceeds to 
construct a homogenized continuum description, then the elastic displacement fields (given in 
dU)) and the energy of stepped surfaces (given in (|l6|) ) are recovered. This procedure is very 
instructive, as it will allow us to give physical interpretation for parameters such as (3i and 
in the continuum description in terms of quantities that determine the elastic fields of discrete 
atomic steps. As we will see, these quantities depend on bonding of the atoms near the step 
edges and therefore have to be determined using an atomic level calculation or, alternatively, 
from experiments. In what follows, we will first look at the elastic displacement fields and the 
corresponding strain energy stored in the crystal with a traction-free stepped surface. As noted 
earlier, the contribution of the strain energy to the total free energy enters through the non- 



local contribution in (16). We will then proceed to analyze the local part of the free energy that 



depends on the local slope of the surface. 
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A.l Non-local parts of the free energy and elastic displacement fields 

In their pioneering work on step interactions, Marchenko and Parshin (1980) analyzed the strain 
field of a surface step by viewing the step as a point dipole on the surface of strength D = 
Diei + 1)262, with components in both the coordinate directions perpendicular to the step 
which lies along the e3-direction in the present case. Since the component of the moment 
should vanish in equilibrium, they concluded that D2 = ±TQhs, where Tq is the surface stress and 
+(— ) denotes a positive(negative) step0 respectively. This value of D2 ensures that the moment 
created by the surface stress is compensated. The other component, Di, depends on details of 
atomic bonding around the step edge and has to be obtained from an atomic level calculation 
or from experiments. The 1-component is a dipole without moment, whereas 2-component is 
a dipole with moment, as noted. If we focus attention on crystal structures that are invariant 
under rotations of 180° about the 62 axis, it is clear that a positive step becomes a negative 
step under such rotations. While the 2-component flips its sign under such a rotation, the 1- 
component does not. Therefore, on such high symmetry surfaces (for example, the (100) surfaces 
of cubic materials commonly employed in strained epitaxy), both positive and negative steps 
have identical values for Di. In addition to these components of the force dipole, for a strained 
crystal equilibrium requirements can be invoked to show that a monopole M = hgCukie^i in 
the ei-direction must be included, in addition to the point dipole, to model the stress fields of 
the step associated with the mismatch strain (Tersoff et al. 1995). We can now use the surface 
Green's function to write the elastic displacement field of an isolated step located at the origin 
as 

n.(x) = ±MGa(x) + D,^^ ± D,^^ (31) 

for a positive/negative step. 

If we consider the surface shown in Fig. ^, with a continuous distribution of steps, the 
surface forces can be modeled using continuum monopole and dipole densities 



= " * and d(x) = ^ 

/lo ax /lo 



dh 



dx 



Do dh 

+ -HIT.'' 



respectively. The displacement fields in the crystal due to of distribution of steps can be written 
as a superposition of the displacement fields of individual steps, or 



00 



, , ^dh(x') DidGiihc- x'ei 
—Gii{x - X ei^ 



-To- 



hs dx' hs dx 

5Gj2(x — x'ei) dh{x' 



dh{x') 



dx' dx' 



dx' 

dx. (33) 



^For a positive step located at xo, the surface height satisfies h{xQ) — h(xQ) = hs, where h{x'^) {h{x^)) is the 
surface height slightly to the right (left) of the step, while for a negative step the relation /i(a;J) — hix^) = —hs 
is satisfied. 
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Once the second and third terms are integrated by parts, it is easily seen that (13 3| ) can be written 
as 



" — 



^^Ga(x - X ei)^ + -Ga(x - X ei)^^Sgn 



+ToGi2[x - X ei) ^^^^ 



dx', (34) 



which is in the form given in (|T4|) with the surface traction f precisely as given by ( [l^ ) if we 
identify /?i with Di/hs- It can also be verified that the contribution to the free energy from this 
displacement field is 

Enon-local = "TT / f • U dxi, (35) 
^ J Si 

which is nothing but the last term in (16) which represents the non-local contribution to the 
free energy. 

To summarize, an interesting outcome of the simple homogenization process that we have 
carried out is that even though individual steps are modeled using both monopoles and dipoles, 
the displacement fields are determined solely in terms of a distribution of monopole-like surface 



tractions given in (IS). The physical reason for this can be understood by appealing to the 
analogous homogenization problem of dipoles in electrostatics or magnetostatics. In a polarized 
dielectric with dipole moment density P(x), the homogenized electric fields can be obtained 
by replacing the spatially varying dipole density with an equivalent charge density pix) = 
— V • P(x); this means that a spatial variation of the dipole density physically corresponds 
to charge accumulation. In the case of steps, the force monopole distribution that originated 



from the dipoles in (13) is seen to be nothing but dd{x)/dx, where d{x) is the dipole density 
introduced earlier. It is also clear that a homogenous distribution of dipoles does not give rise 
to any long-range elastic displacement fields; if we consider a vicinal surface shown in Fig. |l|, 
for regions with size comparable to or larger than the step spacing, symmetry arguments can 
be invoked to show that the average or homogenized displacements in such regions vanish. Of 
course, on length scales much smaller than the step spacing there can be local field fiuctuations 
due to individual dipoles. However, at scales larger than the step spacing, these fiuctuations 
average out to zero and the homogenized field is determined solely by the spatial variations in 
the step density. 

A. 2 Local parts of the free energy 

The local part of the free energy has contributions from two sources, one from the self-energy 
of the force dipoles at the steps and the other due to the interaction of these dipoles with 
the epitaxial mismatch stress. In each of these cases, the interactions arise from the local 
displacement fields near the steps as opposed to the homogenized displacement fields calculated 
in (pj). Below, we consider each of these contributions separately. 
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The self-energy of a force dipole distribution can be obtained from the interaction energy 
between force dipoles on a vicinal surface. If we consider the vicinal surface in Fig. the 
interaction energy of one of the dipoles with all the others can be expressed as 



E, 



dipole 



TT 



11 



(xei) 



=1 + ^2^2^ 



^22(2:61) 
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d2 
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(36) 



where d = hs/ t&nO is the spacing between steps and the factor vr^/G comes from evaluating the 



infinite sum ^ 
using (36) as 



00 , 

i=l 



The self-energy of a continuous distribution of dipoles can be calculated 



E. 



self 



G \h'{xi)\ 



G 



h'{xi)f'dxi, 



(37) 



where the local step spacing d{xi) is expressed in terms of the surface slope as d{xi) = 
hs/\h' {xi)\. The reason that the local approximation works for the self-energy is that the inverse 
square decay of interaction between dipoles guarantees that the contribution of steps that are 
far from a given step are much smaller in magnitude than the contributions from those that are 
close by. In distinct contract, the logarithmic nature of the interaction between force monopoles 
(given in (|l3|)) that arise from the spatial variation of the dipoles, requires us to assign equal 
importance to steps that are both near and far. Hence, this contribution must be treated in the 
non-local approximation. 




Figure 6: Transformation of coordinate axes to calculate the interaction energy of 
the force dipoles and mismatch strain. The coordinate axes in the transformed 
system are chosen so that the I'-axis lies along the vicinal surface and the 2'-axis 
coincides with the normal to the vicinal surface. 

The coupling of the strain field of a step with the mismatch strain gives rise to another 
local contribution to the free energy, which can also be evaluated by considering the vicinal 
surface in Fig. The interaction energy between a positive step and the mismatch strain 
is Diuy ^i' + |-D2|«2',i', where ui'^i' and ui'^2' are displacement gradients in the transformed 
coordinate system shown in Fig. ^. The coordinate axes in this transformed system are chosen 
so that the I'-axis lies along the vicinal surface and the 2'-axis coincides with the normal to the 
vicinal surface. Using the standard rules for transformation of strains, we have 
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U2',i' = 4i + (4-e?i)^ + 0(e'). (38) 

For a continuous surface, the total contribution to the local part of the free energy can be 
obtained by summing the interaction energies of each of the steps on the surface with the 
mismatch strain, or 

Edipoie-strain = I k^^^^(e?i + 2e\^h' {xi)) + |i^2|^^(n° i + (e°2 - el^)h'{xi))\ dxi. 
Js I hs hs ' 

(39) 

Making use of the fact that Di = Pihs and ID2I = Tohg, it can be verified that the terms of 0{h') 
and 0{h''^) in ( p^ ) add up to the local contribution evaluated in ( |39| ) save the factor proportional 
to ^2 1- We do not have an explanation for this difference. 

Within the continuum description the parameter (3i was introduced to include the effect 
of surface stretch on the step formation energy. Comparing the continuum description with the 
discrete step models, we conclude that this parameter is closely related to the intrinsic dipole 
at the step. If there are significant modifications in atomic bonding at the step edges, the 
intrinsic dipole will be large in magnitude. For such steps, strain can significantly reduce their 
formation energy. Atomistic calculations on reconstructed Si steps (Xie et al. 1994, Roland 
1995) have shown large reduction in step formation energies due to strain. The results, however, 
were neither interpreted in terms of the intrinsic dipole nor used in a continuum framework to 
describing evolution of surfaces. Our continuum formalism accounts this crucial effect in the 
description of evolution of strained material surfaces. 

B Energetics of strained islands 
B.l Continuum description 

We will now apply the continuum model to calculate the energy of strained islands. The island 
energies will also be obtained directly by summing up the interactions of individual steps. It will 
be shown that the continuum model not only provides an easy way of calculating the energetics 
but also provides a direct of identifying the contributions to the total energy arising from edges 
and corners. As we will see, extracting these quantities from the discrete step picture is much 
more involved and tedious. 

Let us consider the island made up of N positive and negative steps with base width 2L 
as shown in Fig. ^; the angle made by the facet is given by ^ ~ Nhs/L and the two-dimensional 
volume of the island is A = LNhg- For the sake of simplicity, the surface Green's function of 
an isotropic material will be employed in evaluating the non-local contribution to the energy of 
the island. 

Within the continuum description, the local contribution to energy of the island is 
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2L 

Figure 7: Two dimensional epitaxial island made up of N positive and negative steps 
and with side walls oriented angle 9. The various components of surface tractions 
that are needed to calculate the stress fields in the island and substrate are shown 
in the figure. 



while the non-local contributions to the energy can be calculated from the surface tractions given 
in (jl^). Since the slope on each of the side- walls of the island is a constant, the components of 
the surface tractions are given by 



/i 

/2 



M 



e [e(x -L) + G(x + L)- 2G(x)] + (310 [5{x - L) - 5{x + L)] 



TQe[-25{x) + 5{x -L) + 5{x + L)] 



(41) 



where 0( ) is the unit step function and 5{ ) is the Dirac delta function. The first term in (41)i 
comes from the force monopoles at the steps, while the other terms that act at the apex and the 
edges of the island, as shown in Fig. come from the force dipoles at the step. The different 
components of the non-local part of the island energy can be written as 
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aj5l + 3aro^l 6*^ log(2L//i,), 



(42) 



where the subscripts M and respectively, indicate terms arising from the monopoles and 



dipoles at the steps, and a 

The relative importance of the edge contribution to the island energy, given by (42)2 plus 
(42)3, can be determined by comparing its magnitude with the local contribution given by (^). 
It can easily be verified that the ratios Emd /El and Edd/Ei are equal to 9 and log{A)/A 
within factors of order unity, respectively, where A = A /hi- Since our focus is on shallow 
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islands <C 1) whose dimensions are large compared to atomic dimensions, it is clear that the 
edge contribution can be safely ignored and the energy of the island becomes 



aAe{Ciuje%f\ogA. 



(43) 



If the strain dependence of the first term is ignored, there is an energy barrier to nucleation 
of epitaxial islands (Tersoff and Le Goues 1994). However if the condition /?i + /?ieii < is 
met, for side-wall angles in the range {) < 9 < 9* (see Fig. |3| ) the nucleation barrier is absent. 
Indeed, experimental observations (Sutter and Lagally 2000, Tromp et al. 2000) confirm that 
such shallow stepped mounds islands form as a natural instability during strained heteroepitaxy. 



B.2 Discrete description 

The energy of the strained island in Fig. |^ can be also obtained considering the interaction 
energy of discrete steps that make up the sidewalls of the islands. Using italicized symbols to 
denote energies, the local part of the energy that includes that step formation energy and the 
interaction of the force dipole with the mismatch strain can be expressed as 



8l = 2N [(/3i+/3ieO,)+ro(e02 



(44) 



If we compare this expression with the local contribution in the continuum description, we find 
that the self-energy of the force dipoles does not appear in the discrete framework. The reason 
is that within the discrete model, both the self-energy and the non-local part due the spatial 
variation of the step density have to be evaluated from the dipole-dipole interactions between 
the steps that constitute the island. The total interaction energy, which is a sum of interactions 
between pairs of steps can be written as 
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(45) 

where the steps are labeled with indices i and j and Below, we will evaluate the 

total energy for the island in Fig. |^ by using the elastic Green's function of the isotropic solid. 

The double sum in ( |45| ) can be simplified by rearranging the terms in ( pSD (Shenoy et 
al. 1998, 2000), so that sum over one of the indices can be conveniently evaluated; the step 
interaction energies are then given by 
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When N becomes large, the discrete sums in (^6|)i and (^)2 can be evaluated by converting to 



integrals. The results so obtained are identical to the continuum results in (42). The dipole- 
dipole term given by (^6|)3, in the limit of large N, yields 

£dd = ^^0! + To'] - [ap! + Saro^] 6^ log(2L//i,). (47) 

The first term is identical to the self-energy term of the dipoles in the continuum description, 
while the next two terms agree with the non-local contribution given in (^)3. Thus, the ex- 
pressions for energy of the island computed using both the discrete and continuum descriptions 
are in agreement. 

The above exercise shows that the continuum model provides a direct route to computing 
the energy of stepped surfaces without recourse to tedious evaluation of discrete sums. It also 
provides a direct way of looking at the energies associated with corners and edges. It is difficult to 
obtain these contributions form the discrete approach because the self-energy and the long-range 
part of the dipole interactions are not clearly separated. 

C Convergence of the variational formulation 

The variational formulation is based upon the expression of the surface shape in terms of a 
Fourier series with a finite number of expansion coefficients; refer to (p6t). In this appendix we 
investigate the convergence of the expansion, focusing on the growth a single Fourier mode. The 
goal of this exercise is to show that the amplitudes of the modes with small wavelengths do not 
grow without bound. 

In order to perform the stability analysis of short wavelength modes, the highly nonlinear 
nature of the evolution equations behooves us to consider not just the leading order contribution 
of the energy represented by the negative step formation energy, but also the repulsive step 
interaction energy. It can be shown that the nonlocal contribution, given by the last term in 



(If:), becomes significant only when the wavelengths are large. If only the second and third 



terms in (^) are retained, the evolution of the surface shape consisting of a single mode, 

h{x) = Unit) cos{nkx), (48) 

can be written as 

an{t) = bi + /?ie?i + 2/33a2 fc^n^ SgnK]. (49) 

TT L J 

If j3i + ^leli < 0, the growth of this mode is determined by a competition between the negative 
step formation energy, which promotes the growth of small wavelength modes, and the repulsive 
step interactions. It is evident that a„ achieves the stable amplitude 



iA+M), (50, 
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Because |a„| vanishes as n ^ oo, the use of a finite number of terms in the expansion of surface 
shape in ( p6[) is reasonable. This expectation is borne out by numerical experiments in the full 
problem. 

In the above analysis, we did not include the effects of the substrate. Since the mismatch 
strain is present only in the deposited film, the step formation energy in the substrate is always 
positive. This implies that the amplitude obtained in (|^) represents an upper bound of sorts; 
if the effect of the substrate is included, the stable amplitude would be smaller. 
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